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Abstract 

Introduction: In many practical aspects of ophthalmology, it is necessary to assess 
the severity degree of glaucoma in cases where, for various reasons, it is impossible 
to perform a visual field test - static perimetry. These are cases in which the visual 
field test result is not reliable, e.g. advanced AMD (Age-related Macular Degeneration). 
In these cases, there is a need to determine the severity of glaucoma, mainly on the 
basis of optic nerve head (ONH) and retinal nerve fibre layer (RNFL) structure. OCT is 
one of the diagnostic methods capable of analysing changes in both, ONH and RNFL 
in glaucoma. 

Material and method: OCT images of the eye fundus of 55 patients (110 eyes) were 
obtained from the SOCT Copernicus (Optopol Tech. SA, Zawiercie, Poland). The authors 
proposed a new method for automatic determination of the RNFL (retinal nerve fibre 
layer) and other parameters using: mathematical morphology and profiled segmentation 
based on morphometric information of the eye fundus. A quantitative ratio of the 
quality of the optic disk and RNFL - BGA (biomorphological glaucoma advancement) 
was also proposed. The obtained results were compared with the results obtained from 
a static perimeter. 

Results: Correlations between the known parameters of the optic disk as well as those 
suggested by the authors and the results obtained from static perimetry were calculated. 
The result of correlation with the static perimetry was 0.78 for the existing methods of 
image analysis and 0.86 for the proposed method. Practical usefulness of the proposed 
ratio BGA and the impact of the three most important features on the result were 
assessed. The following results of correlation for the three proposed classes were 
obtained: cup/disk diameter 0.84, disk diameter 0.97 and the RNFL 1.0. Thus, analysis of 
the supposed visual field result in the case of glaucoma is possible based only on OCT 
images of the eye fundus. 

Conclusions: The calculations and analyses performed with the proposed algorithm and 
BGA ratio confirm that it is possible to calculate supposed mean defect (MD) of the 
visual field test based on OCT images of the eye fundus. 

Keywords: Image processing, Measurement automation, OCT, Segmentation, BGA, 
Biomorphological glaucoma advancement, Glaucoma 
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Introduction 

Analysis of OCT images of the eye fundus is a technique known for many years which 
enables to acquire large amounts of information useful for medical diagnosis. Today, 
there are many known profiled algorithms for the analysis of both successive layers of 
the eye fundus [1-20] as well as the choroid layer [21]. Depending on the approach and 
practical implementation, these algorithms operate in the area of texture analysis [1,2], 
algorithms based on segmentation with Canny method [4], SVM [5], mathematical 
morphology [6], active contour [7,8], fuzzy algorithms [9], Markov model [10], robust 
segmentation [11], random contour analysis [12] and others [13-20]. Implementation of 
these algorithms most commonly occurs in the C, C# programming environment if 
they are used in commercially available applications, or, for example, in Matlab if they 
are at the stage of testing and modification [6]. The analysis time of a single OCT 
image of the fundus is several tens of milliseconds to a few seconds with the most com- 
putationally intensive applications [1,17]. In practice, the biggest difficulty is large inter- 
individual variability and a possible high degree of retinal pathology. This translates 
into problems with the designation of successive layers of the eye fundus (RNFL-retinal 
nerve fibre layer, RPE-retinal pigment epithelium, and others) in the following cases 
(Figure 1): 

• The layer (u h u 3 ) is not continuous - Figure la,d, 

• The layer (u 2 , u 3 ) merges with another layer (u 2 ) - Figure lb, 

• One layer ends (u 1 or u 2 ) and another one starts (u 3 ) in the same place - Figure lc. 
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Figure 1 Diagrams showing problems during recognition of layers in an OCT image of the eye 
fundus. The following issues are shown in sequence: a) the problem with the gap between u 1 and u 3 and 
the impact of the additional new layer u 2l b) the problem with combining layers u 2 and u 3 with the layer 
u 1 for 3 pixels (1,4), (1,5) and (1,6), c) the difficulty in deciding which layer u 1 or u 2 should be combined 
with the layer u 3l d) the difficulty in deciding whether the layer u 2 is the missing piece connecting layers u 1 
and u 3 . 
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The mentioned problems can be handled with in a limited way using artificial 
intelligence (knowledge gained from other images) or anthropometric evidence. In any 
case, image analysis and the proposed approach to solving such cases are difficult. 

From a medical viewpoint, the RNFL thickness is essential in the case of the optic 
disk analysis. On the basis of its average thickness and further measurements, other 
parameters such as a cup/disk area ratio or a cup/disk vertical diameter are obtained. 
These parameters have an influence on the diagnosis of an ophthalmologist as well as 
the applied methods of treatment. 

In the case of glaucoma, it is essential to test the visual field. Its measurement is usu- 
ally carried out using a static perimeter. A reliable parameter obtained from perimetric 
data is the average depth of a defect. In many cases, for example, in a number of differ- 
ent additional diseases, it is impossible to obtain this type of measurement or the 
results obtained are not reliable. Then, it is necessary to search for other methods 
which will enable to determine the supposed visual field results. One such method is 
the analysis of OCT images of the eye fundus proposed by the authors. It allows not 
only to determine the severity degree of glaucoma but also to determine the supposed 
MD value of visual field test. It is further shown which selected parameters in OCT 
images correlate, to the greatest extent, with the results obtained from perimetry. 
Moreover, the new algorithm proposed by the authors and the quantitative ratio of the 
optic disk condition are described. 

Material 

OCT images of the eye fundus of 55 glaucoma patients (110 eyes) were obtained 
from the SOCT Copernicus (Optopol Tech. SA, Zawiercie, Poland). Patients 
ranged in age from 37 to 88 years. The obtained images (B-scans) had a resolution of 
MxN = 1000 x 600 pixels (where M - number of rows, N - number of columns) at a 
colour resolution of 8 bits/pixel. For each patient I = 99 B-scans were obtained allowing 
for full 3D reconstruction of the eye fundus image. At the same time, 110 results were 
acquired from the static perimeter OCTOPUS 300Series V 6.07c. All data were analysed 
in a source format. The presented analysis was carried out retrospectively on the basis of 
data obtained during typical medical examinations. The data were anonymised and the 
patients were examined in accordance with the Declaration of Helsinki. 

Methods 

The method of analysis and processing of OCT images proposed by the authors was 
preceded by the analysis of the results obtained from conventional methods. This 
analysis indicates features which should be taken into account when constructing the 
target image analysis algorithm. 

Selection of the features of the image 

As stated in the previous section, a visual field test was performed on 55 glaucoma 
patients. The features of the eye fundus were measured using a Copernicus tomograph. 
The results obtained in the software available for these devices are shown in Figure 2 
and Figure 3. The features obtained on this basis are presented in Table 1. Examples of 
values for each feature for the first 10 patients are shown in Table 2. Each of these 
characteristics will be further, for simplicity, denoted by the symbol "w". In this case, 
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Figure 2 The result for one patient obtained from the Copernicus tomograph. The following images 
can be seen in the corners: fundus reconstruction, a deviation from the normal map, the RNFL thickness 
map and the graph RNFL TSNIT. In the central part at the bottom there is a sample B-scan of the optic disk. 
The central part is dedicated to the ratios calculated by the software. 

/ 



the features from w(l) to w(28) were acquired from the tomograph, and w(29) and 
w(30) from the perimeter. The parameters from w(l) to w(30) were further subjected 
to exploratory factor analysis in order to detect their cross correlation and to remove 
those features which are statistically insignificant. For this purpose, the following items 
were designated in succession: correlation of features w, eigenvalues of the correlation 
matrix, the number of classes using the Kaisers and Cattells criteria, rotation of the co- 
ordinate system, correlation of individual features with the newly created classes. The 
results (for the whole group of patients) are shown in Table 3 whereas the correlation 
chart of various features with three classes are shown in Figure 4. The analysis of 
Table 3 and Figure 4 suggests several conclusions useful for further analysis: 

Class 1 significantly (above 0.5) correlates with the features w(16), w(13), w(6), w(5), 
w(8), w(3), w(15) and w(17) and is in a negative correlation with the features w(19), 
w(7) and w(4). It means that the cup/disk diameter and the cup diameter are a 
representative of the class 1. The ratio of the cup vertical/horizontal and the area 
and volume of the rim are in the inverse correlation. 

Class 2 significantly correlates with the features w(22), w(20), w(24), w(25), w(23). 
These features are responsible for the RNFL parameters. 

Class 3 correlates with the features w(10), w(2), w(12), w(ll) and w(14) which 
are responsible for the parameters of the disk diameter in the vertical and hori- 
zontal axes. 

Analysis of this division into classes, their correlation with the results obtained from 
perimetry and interpretation of individual features allow to specify the features of OCT 
images necessary for further analysis. The particularly important features of OCT 
image are: disk diameter and its relation to cup as well as the RNFL parameters in the 
area. This information is the basis for the design of a profiled image analysis algorithm, 
correlation analysis with the results of perimetry and for suggesting a new quantitative 
ratio of the optic disk condition. 
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Figure 3 The result for one patient obtained from the static perimeter. From the top: data on patient 
and study, a colour picture of the quality of vision in the 30°, numerical data of quality of vision in various 
fields expressed in decibels: values, comparisons, corrected comparisons, probability and corrected 
probability. In the lower right corner there are two most important parameters, namely MD (Mean Defect) 
and LV (Loss Variance). The analysis did not take into account LV value, because there was no significant 
correlation between LV and any of the OCT parameters. 



Pre-processing 

Pre-processing of images is related to automatic unpacking of the source file with 
the extension *.oct and thus acquiring images in the *.jpg format and calibration 
data. These data mainly include information recorded by the tomograph in the 
info.ini file, for example, the position of subsequent B-scans in space, the numbers 
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Table 1 Image features obtainable with the known software of the tomograph and the 
static perimeter 
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of A-scans and B-scans, date of test, the patient data and others. The sequence 
of images L GRAY {m,n,i) obtained on this basis is further subjected to filtration 
with a median filter whose mask /zj is sized M^i x N^i x hi =3x3x3 pixels. The 
mask size was chosen arbitrarily taking into account the size of distortions and 
artefacts entering the optical path. The filtered image L M {m } n } i) (Figure 5) is sub- 
jected to further preliminary transformations. These include determination of the 
RNFL. 

Processing - RPE, ILM, RNFL, cup, disk 

In order to determine the RNFL, the RPE layer (retinal pigment epithelium) is 
pre-determined automatically. The RPE layer is the most characteristic and the 
brightest layer in an OCT image, therefore, it is used to determine the target 
RNFL. Additionally, the RPE is used to determine the disk boundary. For this purpose, 
every n-th column of the image L M) each A-scan, is analysed. Analogously to the method 
described in [21,22], the position of the maximum brightness for each column is pre- 
determined, i.e.: 

L RPE (n, i) = arg[ max (L M (m, n, i)) ) (1) 

\me(l,M) / 

where: 

m,n - coordinates of the rows and columns of the matrix me(l,M) and wg(1,A0- 
The formula (1) can be directly applied only when for all the analysed rows of each 
column there is only one maximum brightness value. As shown in [6], in practice it 
occurs in about 80% of the cases analysed at a resolution of 8 bits per pixel. On this 
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Table 2 Sample results for the first 10 patients obtained from the tomograph and 
perimeter 
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basis, two essential elements are calculated: disk and ILM layer. The disk location, the 
image L disk , is here calculated as the area: 



( 1 if L RPE (n, i) > p m • mezn{L RPE {n, i)) 
Ldisk(n,i) = < nj (2) 

0 other 



where: mean - mean value, 
p m - threshold equal 2/3. 

The value of 2/3 (p m ) was chosen arbitrarily and is associated with the max- 
imum acceptable shift of the RPE layer in the area of the optic disk. The other 
element, namely the ILM layer area, was calculated as the area from the first pixel 
of each row to L RPE (n,i). It is very easy to perform binarization in this area using 
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Table 3 The results (for the whole group of patients) obtained from factor analysis for 
the three classes (bold >0.5, underlined >0.7) 
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Otsu method [23], thus setting the binarization threshold p r . The resulting binary 
image L B (m,n,i) is: 

L (m n i) = i 1 $ iL M (m, n, i) > p r ) a (m < L RPE (n, i)) ^ 
By ' ' ' \0 others ^ } 

Based on the binary image L B (m,n,i), the edge of the ILM layer was defined as: 

oil u{m '"' >) = 1 <«> 

L ILM (n,i)= min (L s (m,n,i)) (5) 

we(l,M) 

The waveform of L ILM {n,i) and the RPE boundary, i.e. L RPE (n,i), are the basis for fur- 
ther analysis of the RNFL layer and cup area. The area between the layers ILM and 
RPE is further subjected to an operation which involves determination of the edge 
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Figure 4 The result of factor analysis of 28 features and their division into 3 classes. Each of the 28 
features is represented in this plot by a vector, and the direction and length of the vector indicate how 
each feature depends on the underlying factors (class). The orthogonal nature of the majority of the 
features, in particular w(16), w{22) to w(25) and w(1), w{2) w(12) as well as w(1 1), is clearly visible. Individual 
observations are marked in red. On the right side in the corners the characteristic areas have 
been enlarged. 




Figure 5 The result of 3D reconstruction of the image L M [m f n,i) and the image L M {n,i). The first image 
is shown in an artificial colour palette in the Cartesian coordinate system. The characteristic layers RNFL and 
RPE, optic disk, and the distribution of blood vessels are visible. This image is the basis for further analysis 
and verification of practical usefulness of the algorithm proposed by the authors. The other image is the 
brightness of the pixels in the matrix system with /' rows and n columns showing the exact location of the 
vessels and the optic disks. 
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using two orthogonal Prewitt masks h h and h v sized 3x3 pixels. Filtration results in 
the image L G (yn,n,i), i.e.: 

L G (m, n, i) = < \J L Gh(m, n, if + L Gv (m, n, if if (m > L ILM (n, i)) A (m < LRp E {n, i)) 
y 0 other 

(6) 

where: L Gh (yn,n,i) and L Gv (m,n,i) are the result of convolution of the image L M (yn,n,i) 
with the masks h h and /z v . 

The RNFL is characterized by a contour visible in each column of the image L M (yn,n,i) 
as a transition from a light area (below the ILM) to a dark one (the other layers). It is 
therefore the first edge below the ILM for L G (yn,n,i) > 0, i.e. (analogously to (4), (5)): 



t (™ i/i i\ J m l f { L c{m, n, i) > 0) a (vn > L ILM (n, i)) a [m < L RPE (n, i)) 
L R {m,n,i) = ^ M ^ 

(7) 

L RNFL {n,i)= min (L R (m,n,i)) (8) 

wg(1,M) 

The protection (m > L ILM (n,i))A(rn < L RPE (n,i)) is additional and prevents detec- 
tion of the RNFL on the edge of the analysed object with the background. The 
resulting edge has a number of disadvantages. The biggest one is the lack of con- 
tinuity caused by both occurring shadows as well as insufficient contrast in the 
place of its detection. For such cases, the method of active edge correction was 
suggested. This method involves adjustment of the position of the layer L RPE (n,i). 
For each point of the layer L RPE (n,i), the mean brightness value above and below 
L RPE {n,i) in the areas L u and L d is determined - Figure 6. The new corrected pos- 
ition of the layer L RPE \n,i) is determined as the maximum difference between the 
average values in the areas L u (kud,n,i) and L d (Aud,n,i). These areas are shifted in 
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Figure 6 The results of the dynamic adjustment of the contour position. On the left side, the idea of 
the dynamic contour position adjustment has been shown. The difference between the mean values in the 
areas L u and L d provides information about the correct position of the point. The greatest difference 
between the areas L u and L d is searched for in a selected area. The right side shows an example of the 
input waveform of L RNFL {n,i) in red and the result of the position adjustment of 

Lrnfl (^/O ' n green. 



Koprowski et al. BioMedical Engineering Online 2014, 13:16 
http://www.biomedical-engineering-online.eom/content/13/1/16 



Page 12 of 18 



the column axis in the range Aud = ±med(\L ILM (n f i)-L RPE (n f i)\) t The adjusted pos- 
ition of the RNFL is therefore calculated as: 



L W(M) = ar g max(\L M (Aud,n,i)-L d (Aud,n,i)\) (9) 

m \ Aud J 

The resulting new position of the layer L RNFL (n,i) is the basis for further analysis. On 
the basis of the positions of the ILM and RPE layers, i.e.: L ILM (n,i), L RPE (n,i), the cup 
area L cup (n,i) is determined. This area is calculated as the intersection of the boundary 
of L disk (n,i) at the offset height equal to 150 um with the layer L ILM {n,i). 

As evident from the exploratory factor analysis presented in the previous section, 
apart from the RNFL thickness, also limitation of the analysis area - ROI (region of 
interest) is vital. The ROI ranges from 2-r disk to 3-r dlsk in an OCT image and in 
order to determine it, it is necessary to determine the position of coordinates of 
the centre (n d0 ,i dc ), i.e.: 

E;=iE„=iW>m) 

l do = w v-^N ; — — y 11 ) 

E*=lE„=lW>M) 

The value of r disk is the average radius calculated from the area of a circle: 



rja . J ZLELw> (I2) 



The calculation of the centre of gravity is followed by the conversion of the ROI to 
the polar coordinate system, i.e.: 

6 = atan2{n-n d0 , i-i d0 ) (13) 

P = \J (n-ridof + (i-idof (14) 

The function atan2 returns 0 the same size as n and i containing the element-by- 
element, four-quadrant inverse tangent (arctangent) of the real parts of n and i. The 
obtained results, the ROI in the polar coordinate system L RNFL \Q,p), are shown in 
Figure 7. The proposed algorithm is shown globally in a block diagram form in 
Figure 8. On the basis of the presented ROI, the following results (features) were 
obtained: the disk diameter and its relation to the cup and RNFL parameters. 
Their thorough analysis and correlation with the results from perimetry are pre- 
sented in the next section. 



Results 

The obtained result L^pi XOyp), the ROI in the polar coordinate system (Figure 7), is 
the basis for further analysis. Based on the presented exploratory factor analysis and 
the proposed method of image analysis and processing, the features representing 
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TE 

Figure 7 The results of the region of interest (ROI) separation in the polar coordinate presented in 
Cartesian coordinate system. The present region of interest is the basis for further analysis related to the 
determination of the features associated with the area change in the RNFL thickness and the percentage 
of the area cup/disk ratio. These features are the basis for determining their correlation with the results of 
perimetry. The areas TE, SU, NA, IN denoting the following areas respectively: temporal, superior, nasal, 
inferior have been marked on the axis 6. The cup boundaries are marked in red and disk boundaries 
in green. 



different classes (Figure 4, Table 3), associated with L RNFL (Q,p), were formulated 
in the following way: 
w g (l) - the relative difference in the positions of centres of gravity calculated as: 



\j (n d Q-n c o) 2 + (ido-icof 



w g (l)=l (15) 

' disk 

where: 

n do> ido - coordinates of the centre of gravity of the disk, 
n co> ico - coordinates of the centre of gravity of the cup, 

w g (2) - relative, inverse, minimum distance between the boundaries of the cup and 
disk i.e.: 

Wg (2) = l- r f± (16) 

i disk 

where: r min - minimum distance between the boundaries of the cup and disk. 

w g (3) - RNFL thickness deviation in the ROI (L RNFL ROI (6)) for particular areas: 
temporal (TM), superior (SU), nasal (NA), interior (IN) from the standard thick- 
ness changes: 



wJ3) =— V 



Lrmfl(®) LYnfl(®) ^7) 



Lrnfl(®) 
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Figure 8 Block diagram of the algorithm proposed by the authors. The patient undergoes a 
fundus examination with the use of a tomograph and their visual field is measured with a perimeter. 
Then, the significance and correlation of the features obtained are calculated, thereby forming three 
classes. They are the basis for proposing a new analysis algorithm. After automatic determination of 
the layers RPE, ILM, RNFL and the cup and disk, the features w g {\), w g {2) and w g {3) are determined. 
In the final stage, the proposed BGA ratio and correlation with results obtained from perimetry 
are calculated. 



where: L RNFL (6) - the reference value of the RNFL thickness (the graph highlighted 
in green in Figure 2). 

The mentioned features are components of the proposed biomorphological glaucoma 
advancement ratio (BGA) which is defined as: 
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Wbga = I - K(l) + w g (2) + ^(3)) (18) 
After value substitution and simplification: 



Wbga — 




Figure 9 presents a list of various types of cup arrangements relative to the disk and 
the impact of their sizes on the values of the features w g (l), w g (2) and w g (3) as well as 
w BGA . In addition, for comparison, the values of the DDLS (Disk Damage Likelihood 
Scale) were also shown [24]. A change in the position of the centre of gravity of the 
cup relative to the disk is expressed by the feature w g (l). The closer the centres of grav- 
ity of the cup are to each other, the closer the feature w g (l) is to zero. If the centres of 
gravity of the cup and disk are far apart, the extreme position of the cup is on the edge 
of the disk, then the value of w g {l) is close to one. In turn, the feature w g (2) is a meas- 
ure of centricity of the distribution of cup edges relative to the disk. When the edge of 
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Figure 9 The impact of cup arrangement with respect to the disk on the values w g and w BGA . 

A summary of various typical cup arrangements relative to the disk and the effect of their size on the 
features w g 0), w g {2) and w g {3) as well as w BGA and, in addition, the values of the DDLS, have been shown. 
A change in the position of the cup centre of gravity relative to the disk is expressed by the feature w g {]). 
In turn, the feature w g {2) is a measure of the centricity of the cup edge distribution with respect to the disk. 
wJ3) is the mean value of the RNFL thickness deviation from the standard. 
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the cup is substantially spaced from the edge of the disk, the value of w g (2) is close to 
zero. When the edges of the cup and the disk overlap in any area, then the value of 
w g (2) is close to one. The situation is similar in the case of advanced glaucoma where 
r min is close to zero and w g {2) is close to one. A change in the value of w g {3) is related 
to a deviation with respect to the standard waveform of the RNFL thickness (Figure 2). 
The greater this difference is, the larger the value of w g (3) is, and vice versa. The 
features w g (l), w g (2) and w g (3) are summed up and thus determine the ratio w BGA . The 
ratio w BGA reaches values close to one if all three features w g indicate pathology. In the 
case of the right features w g , the ratio w BGA is equal to zero. 

Based on the proposed ratio w BGA) its correlation with the results of perimetry was 
determined. In addition, it was compared with the results obtained only on the basis of 
known software (features w(l) to w(29)). The following results were obtained: 

• Correlation of classes 1, 2, 3 (obtained from features w(l) to w(29)) with the results 
obtained from perimetry is 0.78, 

• Correlation of the proposed ratio w BGA referred to the described method of image 
analysis and processing with the results of perimetry is 0.86 - Figure 10. The 
function describing the dependence of w(29) (denoting the MD value from 
perimetry) on w BGA takes the following form: 

w(29) = 24.wbga-9.34 (20) 

The difference in the correlation results stems from the profiled method of image 
analysis and processing and the measurement methodology of the described issue. 
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Figure 10 Graph of the dependence of w{29) as a function of changes in the ratio w BGA . The graph 
shows in blue the results obtained for individual patients. The result of linear regression is marked in green. 
The function describing the dependence of w(29) on w BGA takes the following form: w(29) = 24-w BGA -934. 
On this basis, when the ratio w BGA is calculated, it is possible to calculate the expected result of perimetry. 
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Summary 

The paper presents a new methodology for image analysis and processing profiled for 
the optic disk analysis with a particular emphasis on its correlation with static perim- 
etry. The calculations and analyses performed with the proposed algorithm confirm the 
possibility of prediction of the visual field result based on OCT images of the eye fun- 
dus. The analysis time of a single OCT image of the eye fundus with the use of the pro- 
posed algorithm on the computer Core i5 CPU M460 @ 2.5GHz 4GB RAM does not 
exceed 0.2 s. Additionally, a new ratio reliable with respect to the quantitative assess- 
ment of the optic disk condition was proposed. On its basis, a comparison of the 
results of the correlation between the known methods and the methodology described 
in this paper was made. 

The present paper does not exhaust this interesting subject. The proposed method of 
image analysis and processing is one of the possible approaches to solving this type of 
task. Certainly, other methods of image analysis and processing can be used here 
[25,26] such as texture analysis [27,28] or other methods applied to solving similar 
issues [29-31]. 

In further studies, the authors intend to test practical usefulness of the proposed 
algorithm and the ratio w BGA . The results obtained for large inter-individual variability 
(diversity of shapes and topology of the optic disk), repeatability of the data and the 
impact of technological diversity and device parameters (tomograph and perimeter) on 
the result are particularly vital here. 
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